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Abstract 

We present the first direct evaluation of AI = 3/2 K — > tt-k matrix elements 
with the aim of determining all the low-energy constants at NLO in the chiral 
expansion. Our numerical investigation demonstrates that it is indeed possible to 
determine the K — ► tttt matrix elements directly for the masses and momenta used 
in the simulation with good precision. In this range however, we find that the 
matrix elements do not satisfy the predictions of NLO chiral perturbation theory. 
For the chiral extrapolation we therefore use a hybrid procedure which combines 
the observed polynomial behaviour in masses and momenta of our lattice results, 
with NLO chiral perturbation theory at lower masses. In this way we find stable 



results for the quenched matrix elements of the electroweak penguin operators 
(/ = 2(vr7r|O 8 |^0)_= (0.68 ±0.09) GeV 3 and i=2^\0 7 \K Q ) = (0.12 ± 0.02) GeV 3 
in the NDR-MS scheme at the scale 2 GeV) , but not for the matrix elements of 
O4 (for which there are too many Low-Energy Constants at NLO for a reliable 
extrapolation). For all three operators we find that the effect of including the NLO 
corrections is significant (typically about 30%). We present a detailed discussion 
of the status of the prospects for the reduction of the systematic uncertainties. 



2 



1 Introduction 



The importance of a theoretical understanding of K — > nn decays is underlined by 
the recent measurements of a non-zero value for the e'/e parameter [1,2] (the first 
confirmed observation of direct CP violation) as well as the long-standing puzzle of 
the AI = 1/2 rule. Within the Standard Model of particle physics, the major diffi- 
culty is to quantify the non-perturbative QCD effects. Lattice simulations are being 
used in attempts to overcome this problem, but the results obtained in this way cur- 
rently have large systematic errors. A major source of uncertainty is induced by the 
chiral extrapolation. Lattice simulations are performed with u- and c?-quark masses typ- 
ically in the range m s /2 < m u ^ < m s (where m s is the strange-quark mass), and the 
physical amplitudes are obtained by combining the lattice measurements with Chiral 
Perturbation Theory (^PT). Previous calculations of K — > titt matrix elements were 
performed at leading-order (LO) in the chiral expansion. Following some early attempts 
to evaluate K — > irir matrix elements directly [3-6], lattice estimates of the AI = 3/2 
K — > irir matrix elements (7T + 7r \O^g\K + ) [7-11] have generally been obtained by cal- 
culating the K — > 7T matrix elements (tt + \Oj^\K + ) and using the soft-pion theorem 
(— i)(ir\0\K) / f n ~ (ttit\0\K) [12]. An exception to this has been the computation of 
the K — > irir matrix element of the operator + (defined in eq. (4) below), with all three 
mesons at rest [13]. 

We have previously proposed improving the precision of lattice determinations of 
the decay amplitudes by computing directly K — > tttt matrix elements at unphysical 
kinematics and combining them with xPT at next-to- leading order (NLO) [14,15]. By 
comparing the behaviour of the computed matrix elements as functions of the mesons' 
masses and momenta with the predictions of chiral perturbation theory at NLO, we can, 
at least in principle, determine the low-energy constants and hence evaluate the physical 
amplitudes. In this paper we perform an exploratory numerical study to investigate the 
extent to which such a programme can be carried out in practice. In particular, we 
evaluate the matrix elements (7r + 7r°|(9j|i^ + ) for the following AI = 3/2 operators 0{. 



O^ 2 = (s a d a ) L (uf3Up) R - (dpd^R + {s a u a ) L {updi3) R 



i^ada) L 



(d p d 



+ \s a u c 



PJL, 



OT = (s a d p ) L 



[upu a ) R - (dpd 



a)R 



{s a up) L {upd a ) R , 



(1) 
(2) 
(3) 



where a and (3 are colour indices and (iPi4>2)l,r = ^17^(1 T 7s)V'2- 



The matrix elements of the AI = 3/2 operator O4, first introduced in ref. [16], give 
dominant contributions to CP-conserving K — > tttt decays. O4 is proportional to the 
AI = 3/2 component of the operator 

0+ = - [{s a d a ) L (u a u a ) L + (s a u a ) L (u a d a ) L ] , (4) 
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whose matrix element, (tt + it \O + \K + ) = (7r + 7r°| C>4| K + )/3, is an important theoretical 
input to the study of the AI = 1/2 rule. A consequence of chiral symmetry is that 
this matrix element can be related to (K°\(sd) l{scI) l\K°) [17], which contains the non- 
perturbative QCD effects in K°—K° mixing (i.e. the two matrix elements are given in 
terms of the same low-energy constants). 

O^ 2 and O^ 2 are proportional to the AI = 3/2 component of the electroweak penguin 
operators, 

3 3 
7 = -(s a d a ) L J2e q (qpq/3)R and O s = -{s a d p ) L ^ e q {q p q a ) R . (5) 

1 q 1 q 

These operators, although suppressed by a factor of a cm /a s compared to QCD penguin 
operators, can contribute significantly in the theoretical prediction of e'/e [18,19]. Recent 
estimates [20-24] suggest that Og has a significant negative contribution to e'/e. 

In principle, by using finite-volume techniques, it is possible to calculate K tttt 
amplitudes directly from lattice QCD [25,26], without using chiral perturbation the- 
ory. However, in order to obtain the physical amplitude from simulations with periodic 
boundary conditions for the fields, one requires lattice sizes of order 6 fm or larger, which 
is numerically very demanding 1 . In addition for AI = 1/2 decays, the lack of unitarity 
and enhanced long-distance effects in quenched and partially quenched QCD suggest that 
one needs to use full QCD to determine these matrix elements [29-31]. This problem is 
much less severe for the AI = 3/2 decays studied in this paper. For matrix elements 
as complicated as those contributing to K — > nn decay amplitudes, unquenched simula- 
tions on sufficiently large volumes and on fine lattices are unlikely in the near future. It 
will therefore continue to be necessary to combine with lattice simulations for the 
foreseeable future. 

In this paper we present the first lattice study of the matrix elements (7r + 7r°|(9j|A' + ) 
(Oi = 04, g 2 ) at next-to-leading order (NLO) in the chiral expansion. We follow the 
strategy proposed in refs. [14,15]. Under the S77(3)l x S77(3)r flavour transformation, 
the operator C 4 is in the (27l, Ir) irreducible representation and O^g (and hence 0?g) 
are in the (8l, 8r) irreducible representation. This means that the chiral expansion for 
the matrix elements of O4 starts at 0(p 2 ), while that for (9 7 8 and 0^ 8 2 starts at O(p ). 
In sec. 2 we recall the demonstration that by allowing one of the pions to carry non-zero 
spatial momentum we are able, in principle at least, to determine these matrix elements 
at NLO in the chiral expansion [14,15]. 

Although we were able to determine the matrix elements (7r + 7r°|(!?4|A' + ) at the quark 
masses at which the simulations were performed, we were unable to perform the chiral 
extrapolation with any degree of confidence. For this matrix element there is one low- 
energy constant at lowest order and six more at NLO and it was not possible to determine 

1 However, see refs. [27] and [28] for the possibility of reducing the required size of the volume by 
using different boundary conditions. 
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so many parameters with sufficient precision by using our lattice data. Nevertheless, our 
study already reveals some qualitative feature of this matrix element, which will be 
discussed in detail in sec. 5. For (7r + 7r° \0^\K + ) there are fewer low-energy constants 
and we were able to estimate the physical matrix elements. Our best estimates for these 
matrix elements renormalized in the NDR-MS scheme at 2 GeV are (see sec. 7), 

I=2 (nn\0 8 \KXt° s = 0-68 (6) (4) (5) GeV 3 

rM™\0 7 \KX™ s = 0-12 (1) (1) (1) GeV 3 . (6) 

The first error is statistical, the second is due to the uncertainties in the non-perturbative 
matching to the continuum renormalization scheme and the third is due to the chiral 
extrapolations. In addition to the errors quoted above, there are uncertainties due to 
quenching, possible scaling violations and finite-volume effects which are more difficult 
to estimate reliably. We discuss these in sec. 6. 

The plan for the remainder of the paper is as follows. In the next section we explain 
our strategy for the determination of all the low-energy constants which are required 
for the evaluation of the physical amplitudes at NLO in the chiral expansion. The 
parameters and details of our simulation are presented in sec. 3. In sec. 4 we explain 
our procedure for the non-perturbative renormalization of the lattice operators. The 
details of our analysis and a discussion of the systematic uncertainties are presented in 
sees. 5 and 6 respectively. Section 7 contains the main numerical results of this study 
and a comparison with those from previous calculations. Preliminary reports on this 
work have been presented in refs. [32,33]. We present our conclusions in sec. 8. There 
is one appendix in which we study the chiral behaviour of the K — > it matrix elements 
obtained by the CP-PACS collaboration [11] and argue that different choices of ansatz 
for the chiral extrapolation can lead to significantly different results. 



2 Strategy for the Determination of the Physical 
Amplitudes 

In this study we follow the strategy presented in refs. [14,15] for the determination of the 
physical amplitudes from lattice simulations combined with %PT at NLO. In order to be 
able to extract all the low-energy constants at NLO, we must evaluate K — > tttt matrix 
elements for a range of unphysical masses and momenta. We choose to use the SPQR 
kinematics, in which one of the final-state pions has a non-zero spatial momentum and 
the other two mesons are at rest. 

The chiral expansion for the matrix elements A^ S pq R = (7r + 7r°|04|A' + ) and A^spqr = 



3 



ir + 7r°\Oj^\K + ) in the SPQR kinematics, to NLO and in infinite volume is: 



A4 (4) 



6 /2 f 1 

a(E w M n + -M K (K + M n ) + [chiral logs]$ QR 



+4&M* + (4& + Wt)E«MI + (ft - ft + f3 7 )M*M K 

+(ft - ft + ft + 2p 22 )E 7T M*M K + (-4ft + 8ft 4 )M 2 M 2 + (ft + 2ft)£ 7r M| 
+(-2/3 5 + 4ft + Af3 22 )E n M n M 2 K + (ft + 2{3 7 )M n M 3 K 



+(-l6(3 2i )E 2 M 2 + 2(3 22 E 2 M n M K + 8(3 2i E 2 M 2 K 



QV2( 



a{ [f7 ff M ff + -M x (E n + M 7r )J [chiral logs] /8 (M w , Af*) + [chiral logs]^ QR } 



24L 4 16L 5 \ „ „ ,„ 



+4ftM* + 4ft + 2ft 



/ 2 / s 



+ 



+ 



+ 



ft-ft + & 



/? 4 - /5 5 + /?7 + 2/? 22 - 



— - y) * 

12Z^4 8-L5 



/ 2 / 2 



E V M$M K + (-4ft + 8(3 M )M 2 M 2 K 



^ + 2^ 7 -^-^)^M| 
-2/? 5 + 4/3 7 + 4/3 22 - 



/ 2 / 2 
_ ^4 _ 41^ M ^ M 3 



48/^4 8Z/5 \ „ , , „ ,n 



/ 2 f 



litivj. k 



l6(3 2i )E 2 M 2 + 2(3 22 E 2 M„M K + 8(3 M E 2 M£ 



(7) 



(7,8) 
SPQR 



4^2 

/*/ 2 



7 



(7,8) 



(1 + [chiral logs]f P w Q p R ) + 2(5^ 8) + 8™ + 25^) M] 



+ 2 (2(4 7 ' 8) + ^ 7 ' 8) ) + 4 7 ' 8) K - (5f< 8) + 8^)m^ 



| 7 (7 ' 8) (l + [chiral logs]|WP R + [chiral logs] /S (M ffJ M K 

+2(e 8) + # 8) + 24 7 ' 8) - ^ - jf) Ml 
+\(8™ -8™ -25™){M^E„)M K 



2 

+2(2(6™ + $' 



4 ? ' 8) ~ l ^-^)Ml- (5^ + ^ 8, )M W ^} , (8) 



where [chiral logs]gpq R and [chiral logsjf^q^ are calculated in [15] in both full QCD 
and quenched QCD (qQCD). With the exception of /, which is the pseudoscalar decay 
constant in the chiral limit 2 , the quantities appearing in eqs. (7) and (8) all correspond 
to the masses and energies used in the simulations (M w is the pion mass, Mk the kaon 
mass and the energy of the pion with non-zero spatial momentum). All dimensionful 
quantities are taken to be in lattice units, a, 7( 7,8 ) and #j are the unknown low- 
energy constants (LECs) in the chiral expansion for these weak amplitudes. are the 
standard Gasser-Leutwyler constants in the strong chiral Lagrangian [34]. [chiral logs]/3 
represents the chiral logarithms appearing in the renormalization of the factor 1/ fxfl 
and are known for both QCD and qQCD [34,35]. 

The chiral expansion for the physical matrix elements, M,^ y f (that is (7r + 7r°|(94|fr + ) 

and (7r + 7r° | O^'g 2 with the four-momenta of the initial and final states equal), to NLO 
in infinite volume is: 

= - M « + [ Chiral l0g *t 

+(/? 4 -& + 4/3 7 + 2/3 22 )il4 



+(4/3 2 -4/? 4 -2/5 7 -16/? 24 )M^ 

+(-4& + 3/3 4 + & - Wi ~ W22 + IQ(3 2 a)M 2 k MI 



.^5 ja((M£ - M*)[chiral logB] /S (M ffJ M K ) + [chiral logs]^ 



+ ( Pa ~ Ps + 4/3 7 + 2/3 22 - - -± ) M 4 K 



48L 4 8L B 



2 We work in the convention for which the physical value of / ff is 132 MeV. For the physical pion and 
kaon mass we have used the values M v = 137 MeV and Mk = 495 MeV. 
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+ Ufk - 4/?4 - 2(3 7 - 16/3 24 + + yp) 

+ f-4ft + 3/3 4 + & - 2/3 7 - 2/5 22 + 16/3 24 + ^ - ^ J M|m4 , (9) 



M™ = (1 + [chiral logs]™- 

+( - 4 7,8) - 4 7,8) + 2(4 7,8) + 4 7,8) ) + 44 7>8) ) mi 

+ (5f' 8) + 4 7 ' 8) + 4(5f' 8) + 4 7 ' 8) ) + 26^)M^ 

| 7 (7 ' 8) (l + [chiral logs]™ p + [chiral \ogs} P (M n , M 1 

,(7,8) _ ,(7,8) , ofx (Jfi) , ,(7,8)x , ,,(7,8) _ 48^ 81^ 2 



+ 4 7,8) + 4(4 7 ' 8) + 4 7 ' 8) ) + 24 7,8) - ^ - ^H 2 } . (io) 

We use the formulae for [chiral logs] p k ys and [chiral logsj^f for both full QCD and 
quenched QCD (qQCD) presented in [15]. 

In eqs. (7) -(10), we have written the amplitudes using two expressions which are 
equivalent at NLO in the chiral expansion. The difference corresponds to the two ways 
that we incorporate the chiral corrections to l// 3 at NLO. This leads to two procedures 
for obtaining the LECs by fitting -Mspqr to a function of the masses and energies: 

1. we multiply the amplitudes given by the first expressions in each of eqs. (7) and 
(8) by fxfl computed at the same values of the quark masses, and obtain the 
parameters a, the /3's, 7 and the <5's by fitting; 

2. we determine the ratios «// 3 , the /3// 3 's, 7// 3 and the <5// 3 's by fitting the second 
expressions in eqs. (7) and (8). 

Notice that in the second approach, we do not need to determine the Gasser-Leutwyler 

constants L4 and L 5 in order to obtain the matrix elements, because these constants can 

(7 s) 

be absorbed into the LECs and 5; ' via the re- definition: 

At — > Pi = Pi + ~JT'> 
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(4 7>8) + C ,8) ) - (4 7,8) + # 8) y = (4 7,8) + 4 7,8) ) + 



4L 5 



/2 > 



5 t 



■(7,8) 
6 





12L 4 



(11) 



P 



The difference between these two procedures is at next-to-next-to leading order, hence 
beyond the precision we are trying to reach in this work. We find however, that perform- 
ing the fits using the second procedure leads to more stable fits and smaller statistical 
errors. We therefore quote the results obtained using the second procedure. 



3 Details of the Simulation 

The results presented in this paper were obtained from two sets of data. Set 1 with 340 
quenched gauge configurations and set 2 with 480 configurations, both on a 24 3 x 64 
lattice at /3 — 6.0 (which corresponds to an inverse lattice spacing a" 1 = 1.98(6) GeV 
as determined from the masses of the K and K* mesons) but on which different valence 
quark masses have been simulated. The gauge action is the standard Wilson plaquette 
action, and the fermions are described by the Sheikholeslami-Wohlert (SW) action [36] 
with the coefficient csw tuned so that the action does not contain 0{a) artefacts to all 
orders in a s [37] (at our value of /?, csw = 1-769). The 0(a) improvement of the lattice 
four-quark operators requires mixing with operators of higher dimension and is beyond 
the scope of this paper. 

In tab. 1 we list the masses of the pseudoscalar mesons used in our calculations. For 
the pions, for each data set we take the three combinations with a degenerate quark- 
antiquark pair. For the kaons we take all quark-antiquark combinations. For each 
pion the corresponding kaon has one quark or antiquark with the same mass (the kaons 
corresponding to each pion are listed in tab. 1). Since we wish to exploit the strategy 
proposed in refs. [14, 15], the kaon and one of the final-state pions are at rest, while 
the spatial momentum of the other pion is either p = or p = p min = (2ir/L, 0, 0) (all 
physical quantities are averaged over the three directions). 

In order to determine the matrix elements of the lattice operators Oi = O4, 7 
and O^ 2 (where the bar denotes the bare lattice operator) we evaluate the correlation 
functions 



C$\t,p) = (0|P w+ (t fix ,0)P^o(t,p)(5 J (0)P| + (-t fix ,0)|0) , 



# 


«1 


K 2 


aM P 


aE n(P m in) 


meson 


Set 1 (340 configurations) 


1 a 
Id 


D 1 'MzLQ 

U . 1 O^t^tcJ 


D 1 1AAQ 


D 9^71 (\ F\\ 

U. ZJ ( 1 1 1U J 




JJ1U11, KdUIl 


za 


u. loQ^iy 


U. loZDU 


u.oyuz^io ) 




kaon 


3 a 

OCX 


1 3449 


12997 


5177f12"l 




kaon 


lb 


0.13376 


0.13376 


0.3593(13) 


0.4431(16) 


pion, kaon 


2b 


0.13376 


0.13146 


0.4823(10) 




kaon 


3b 


0.13376 


0.12565 


0.7215(10) 




kaon 


lc 


0.13300 


0.13300 


0.4432(9) 


0.5127(11) 


pion, kaon 


2c 


0.13300 


0.12994 


0.5853(8) 




kaon 


3c 


0.13300 


0.12612 


0.7343(8) 




kaon 


Set 2 (480 configurations) 


la 


0.1344 


0.1344 


0.2716(15) 


0.3770(23) 


pion, kaon 


2a 


0.1344 


0.1337 


0.3224(14) 




kaon 


3a 


0.1344 


0.1328 


0.3785(12) 




kaon 


4a 


0.1344 


0.1318 


0.4338(11) 




kaon 


5a 


0.1344 


0.1308 


0.4839(11) 




kaon 


6a 


0.1344 


0.1299 


0.5259(11) 




kaon 


7a 


0.1344 


0.1282 


0.5995(11) 




kaon 


lb 


0.1337 


0.1337 


0.3663(12) 


0.4488(15) 


pion, kaon 


2b 


0.1337 


0.1344 


0.3224(14) 




kaon 


3b 


0.1337 


0.1328 


0.4173(10) 




kaon 


4b 


0.1337 


0.1316 


0.4790(9) 




kaon 


5b 


0.1337 


0.1306 


0.5260(8) 




kaon 


6b 


0.1337 


0.1295 


0.5746(8) 




kaon 


lc 


0.1328 


0.1328 


0.4642(8) 


0.5309(10) 


pion, kaon 


2c 


0.1328 


0.1344 


0.3785(12) 




kaon 


3c 


0.1328 


0.1337 


0.4173(10) 




kaon 


4c 


0.1328 


0.1316 


0.5220(7) 




kaon 


5c 


0.1328 


0.1307 


0.5625(7) 




kaon 



Table 1: The k values and the corresponding non-degenerate pseudoscalar masses Mp 
for the two sets of data. The pions correspond to the degenerate cases la, lb, lc while 
the kaons respectively to {la, 2a, . . .}, {16, 26, . . .}, {lc, 2c, . . .}. -En-(p m in) is the energy 
associated with a pion when it carries the spatial momentum p min of one lattice unit. 
All the errors are statistical and obtained with the jacknife procedure. 
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C?(t,p) = (0|P^o(t fix ,0)P^ + (t,p)A(0)Pl + (-t fix ,0)|0) 



(12) 



where 



Pp(t,p) 



J 



d 6 



—ip-x 



(f) P (x, t), 



(13) 



and (pp(x,t) is the local interpolating operator for the pseudoscalar meson P (e.g. 
<f) n +(x) = d(x) , y 5 u(x)). In the correlation functions in eq. (12), the spatial momen- 
tum p is either zero or p min , while tfi x is the time for which the two pions propagate and 
in our study it is fixed to be 12. In a preliminary study performed on 140 gauge con- 
figurations with pion masses covering the range of the present simulation we considered 
the three cases t^ x = 10, tfi x = 12 and t^ x = 14. Observing no discrepancy we concluded 
that tfbc = 12 is an appropriate choice. 

In addition to the four-point functions in eq. (12), in order to eliminate the matrix 
elements of the interpolating operators P n +,o and Pk+ and to compute the pseudoscalar 
decay constant we need to evaluate the two-point correlation functions 



where A is the time component of the local axial- vector current. 



4 Renormalization of the Operators 

To obtain the physical amplitudes, lattice-regularised operators have to be matched 
to some continuum renormalization scheme (normally the MS scheme) in which the 
Wilson coefficients are calculated. If chiral symmetry is preserved by the regularization, 

3/2 3/2 

OT and 0° 8 ' mix between themselves, while renormalizes multiplicatively. In a 
Wilson-like regularization, such as the one adopted in this work, chiral symmetry is 
explicitly broken and the bare lattice operators 04(a), 0^ 2 (a) and Og^ 2 (a) mix with 
other dimension-6 operators with different chiral properties 3 . Using CVS symmetry it 
can be shown that mixing with operators of different chirality only occurs in the parity- 
even sector [38-40]. Therefore the renormalization properties of the parity-odd AI = 3/2 
operators considered here are not affected by the breaking of chiral symmetry due to the 
lattice fermion action. The bare and renormalized operators are therefore related by 



3 In the SU(2) isospin limit, the operators do not mix with those of lower dimension, as this mixing 
can only occur through penguin contractions, corresponding to AI =1/2 transitions. 



C PP (t,p) = (0|P P (t,p)4(0,0)|0>, 

c PA (t,p) = (o|p P (t,p)4(o,o)|o), 



(14) 
(15) 




(16) 
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and 





f Z 77 (fia) Z 78 (fia)\ f0 3 7 /2 (a) 
yZgrifia) Z 88 (iia) J \0 3 8 /2 (a) 



) 



(17) 



where fi is the renormalization scale. The renormalization constants depend, of course, 
on the continuum scheme used to define the renormalized operators. In the following we 
will consider the RI-MOM and the MS schemes for which the Wilson coefficients at the 
next-to-leading order (NLO) have been computed in refs. [41-47]. 

In perturbation theory at NLO, matching of lattice and renormalized operators re- 
quires a one-loop lattice (and continuum) calculation of the matrix elements of four- 
quark operators between quark states. However, lattice perturbation theory frequently 
has large coefficients and we know from many examples where a non-perturbative deter- 
mination of the renormalization constants is possible that one-loop lattice perturbation 
theory may be insufficiently precise. Moreover, it has been shown in refs. [8,9] that the 
large anomalous dimension and the scheme- dependent constants in Z result in significant 
uncertainties in the determination of the matrix elements of 7 § 2 4 . 

Several methods have been proposed in attempts to improve the accuracy in the 
determination of the mixing coefficients. One of these consists in evaluating the renor- 
malization constants computed in perturbation theory using an effective coupling which 
resums tadpole contributions [48,49] and which is expected therefore to reduce higher 
order corrections; we refer to this as Boosted Perturbation Theory (BPT) (there are also 
more elaborate variations of this approach which we will not consider here). However, 
as will be seen in sec. 4.1 below, since the one-loop coefficients in Z, and particularly in 
Z 77 and Z 88 , are so large, it is difficult to be confident that higher order terms would not 
change the results significantly. It is therefore important to determine the renormaliza- 
tion constants non-perturbatively. We implement such a non-perturbative method, in 
which these constants are determined from the computation of Green functions between 
quark and gluon external states, as proposed in ref . [50] . This method has already been 
extensively applied to the parity-even component of the AS" = 2 and AI = 3/2 opera- 
tors. A detailed discussion can be found in ref. [51]. In the calculation of the physical 
amplitudes we will use both the perturbative and the non-perturbative renormalization 
constants, although our final results will be those obtained with the non-perturbative 
method. 

4 The MS two-loop anomalous dimensions are also large for the electroweak penguins. One can 
also determine the lattice two-loop anomalous dimensions and they are even larger than their MS 
counterparts. 
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4.1 Perturbative matching 



The renormalization constants Z+ and Z defined in eqs. (16) and (17) can be related to 
the one-loop perturbative corrections of local quark bilinears [38]. The general formulae 
relating the renormalization constants of the bilinear operators to those of the four- 
fermion operators in a generic lattice regularization have been derived in [7]. In this 
subsection we determine the one-loop matching coefficients of the four-quark operators 
O4 and Oyg 2 between the lattice regularization which we are using and the MS-NDR, 
MS-HV ('t Hooft-Veltman) and RI-MOM renormalization schemes. We use the matching 
coefficients for the quark bilinears from refs. [52,53] and combine them with the formulae 
which relate the perturbative corrections to four-quark and bilinear operators and which 
were obtained using Fierz identities and charge conjugation symmetry [7,38] 5 . For the 
RI-MOM scheme we find: 

z ri-mom = i + ^L (_41n(/ia) - 37.517) , (18) 



z ri-mom = I l + SH-2 ln(H " 34.604) £ (6 m(H - 12.672) 



f£ (-2.386) 1 + ^(16 ln(/xa) - 65.460) 



In order to convert this result into the MS-NDR scheme or into the MS-HV scheme we 
give the matching coefficients r+ DR,HV and f NDR,HV [46] defined through 

qNDR,HV = M - a>3 ^ r NDR,Hv\ qRI-MOM ^Q) 



^q3/2^nDR,HV _ ^ _ a s(^) ~NDR,HV j ^q3/2^RI-MOM ^1) 

The values of r+ DR ' HV and f NDR ' HV are 



and 



14 , NnR { | + 1 In 2 -2 - 21n2 

- NDR - 81n2 f NDR =\ 3 3 22 

1 2-21n2 -f + |ln2 1 



^ = -2-81n2 f HV = ( "3 + 3 ln2 - 8 - 2 2 ln2 ] . (23) 

-2-21n2 -f + |ln2 1 



The coefficients in eqs. (18)-(19) are very large. In estimating the numerical values 
of the renormalization constants, we try to minimize the higher-order contributions by 

5 Wc take csw = 1 in the lattice perturbation theory, which is the consistent procedure at one-loop 
level [53]. 
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using a boosted coupling constant denned by: 

S. = ^-, where a?) = |— , (24) 

a^ ** is the bare strong fine-structure constant and V is the average value of the plaquette. 
We perform our simulation at (3 = 6.0 for which V = 0.59370 and a s = 0.13404. This is 
the value of the coupling which we use in our numerical estimates of the renormalization 
constants. For example, at \x = 2.0 GeV we obtain 

Z f~ MOM = 0.600 and i™ = [ °' 631 -° 135 ] . ( 25 ) 
+ ^-0.025 0.302 J 

The numerical estimates of the renormalization constants which we obtain using the 
boosted coupling in eq. (24) are in reasonable agreement with those obtained non- 
perturbatively: compare, for example, the results in eqs. (25) with those in eqs. (33) 
and (34) below. However, the one-loop corrections are so large that we would not have 
trusted the values in eqs. (25) had we not also evaluated the renormalization constants 
non-perturbatively (see section 4.2). On the other hand, consistency with the pertur- 
bative calculation gives us confidence that the large deviations from unity (or the unit 
matrix) found non-perturbatively for the renormalization constants are genuine and not 
simply a consequence of lattice systematic errors. 



4.2 Non-perturbative matching 

The procedure necessary to compute the renormalization constants of A J = 3/2 four- 
quark operators non-perturbatively using the RI/MOM method has been presented in 
detail in ref. [51] following the general philosophy of ref. [50]. Recently, an extensive 
analysis of the renormalization constants of both bilinear and four-quark operators for 
the non-perturbatively (9(a)-improved Wilson action in the quenched approximation has 
been presented in ref. [54]. In that study, four different values of the lattice coupling 
were considered: j3 = 6.0, 6.2, 6.4 and 6.45. For (3 = 6.0, two lattice volumes were used: 
16 3 x 52, with 500 gluon configurations, and 24 3 x 64, with 340 configurations. Since 
the data for the large volume and this work's data set 1 are exactly the same, here we 
will briefly outline the renormalization procedure and refer the reader to ref. [54] for full 
details. 

The non-perturbative determination of renormalization constants with the RI/MOM 
method is based on the numerical evaluation, in momentum space, of correlation func- 
tions of the relevant operators between external quark and gluon states. For the operators 
04, O^ 2 and O^ 2 , we compute the Green functions Gi(p) with i = 4, 7, 8 (for clarity, 
we suppress colour and spinor indices) 

Gi(p)= I II d'xke-^-^ + ^-^^M^^^dmM^Mx,)) , (26) 

J k=l 
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with all four external legs at equal virtuality p 2 = fi 2 . The amputated Green functions, 
Aj, are defined by 

Hp) = S-'ip) ST 1 ® G t {p) S~\p) S'\p) , (27) 
where S is the quark propagator 

S{p)= f d 4 xe- ip - x (^{x)^(x)) . (28) 



The RI/MOM renormalization method consists in imposing that the amputated Green 
functions, Aj(p), computed in the chiral limit, in a fixed gauge and at a given large scale 
p 2 = fi 2 , are equal to their tree- level values. In this study we work in the Landau gauge. 
In practice, we use Dirac projection operators, Pj, to implement the renormalization 
conditions 

Z^ 2 (fia) Z ik {iia) Tr [A fc (p) P 5 \ 2 _ 2 = (29) 

V — A* 

in the chiral limit. Z$ is the quark field renormalization constant, which we calculate 
from the quark propagator as in ref. [51]. The projection operators satisfy the orthogo- 
nality relation Tr [A- ^ Pj] = 5ij where A- ^ is the tree-level amputated Green function 
of the operator CV In this work, we use the projector basis defined in ref. [51]. 

For the RI/MOM method to be applicable, the renormalization scale \x must be much 
larger than Aq C d, in order to be able to use the Wilson coefficients computed in per- 
turbation theory, and much smaller than the inverse lattice spacing, to avoid large dis- 
cretization errors, so that formally 

Aqcd < H < a' 1 . (30) 

The precise range of validity of the RI/MOM method for the operators studied in this 
paper is investigated in the following subsection. 

It is important to notice that the validity of the RI/MOM approach relies on the 
fact that non-perturbative contributions to the Green functions vanish asymptotically 
at large p 2 . A possible difficulty in the implementation of the above procedure comes 
from the coupling of the operators to the Goldstone boson [55-57]. In ref. [54] it is shown 
that for parity violating operators, such as the ones considered here, only a single pion 
pole can be presented in the amputated projected Green function, while in the parity 
conserving case, the appearance of a double pole is also possible. The Goldstone boson 
makes the chiral limit of the amputated projected Green functions singular and hence 
the extraction of the renormalization constants becomes difficult and affected by large 
uncertainties. In order to subtract the Goldstone pole before extrapolating to the chiral 
limit, we follow the strategies described in detail in ref. [54] which are based in part in the 
method suggested in ref. [58]. In ref. [54] it is demonstrated that the contribution of the 
Goldstone pole may indeed affect the extrapolation of some renormalization constants 
to the chiral limit but this problem can be overcome by non-perturbatively subtracting 
the singular term. In this work, we confirm the conclusions of ref. [54] and in the next 
subsection we present the numerical results for the renormalization constants. 
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4.3 Renormalization Group Behaviour of the Renormalization 
Matrix 



The renormalized operators, obtained either perturbatively or with the non-perturbative 
method, should match the renormalization dependence of the Wilson coefficients in order 
to obtain physical amplitudes which are independent of fi. Whereas this is automati- 
cally enforced by a consistent use of the perturbative renormalization constants at the 
next-to-leading order, it is useful to check whether the non-perturbative renormalization 
constants also satisfy the expected renormalization group behaviour. Following ref. [54], 
we thus define the renormalization group invariant (RGI) renormalization constants by 

Z RGI (a) = w?Qjl)Z + (jjl) Z RGI (a) = (w-Y^Z^) , (31) 

where contains the scale dependence of the Wilson coefficients. At NLO we write 



1 + S!M r 

47T 



(32) 



with a similar expression for w^ 1 ^) in terms of J + and 7+ . The values of J + and J 
have been computed at NLO in [47]. As mentioned above, if we compute the quantities 
in eq. (31) by using the non-perturbatively determined renormalization constants then 
the scale independence is not guaranteed. In fig. 1 we show the numerical results for the 
quantities in eq. (31) as functions of \i. 

Although the computed values of the Z RGI in fig. 1 are independent of fi to a rea- 
sonable precision, we see that especially for Z%% the plateau is worse than in the other 
cases. This could be due either to lattice artefacts or to the possibility that higher order 
contributions to the RG evolution may be significant. The detailed study in ref. [54], 
using data from simulations at four values of the lattice spacing, suggests that the lattice 
artefacts are not large and interpret the deviation from scaling as being due to N 2 LO per- 
turbative corrections. As emphasized in ref. [54] it would be desirable to control better 
the scaling behaviour by evaluating the N 2 LO anomalous dimensions and/or performing 
a step-scaling analysis. 

We extract Z R and Z RGI by fitting each of the plateaus in fig. 1 to a constant in 
the interval 2.1 GeV < fi < 2.4 GeV. We then use eq. (31) to obtain the value of the 
renormalization constants at the required scale. A standard scale at which the matrix 
elements are matched with the Wilson coefficients is \i — 2 GeV. At this value of /j,, the 
values of the renormalization constants in the RI-MOM scheme are 

z RI-MOM = 601 (3) (33) 



Iri-mom = I 0-664(2) (11) -0.142(1) { _ 2j 

-0.055(1) (ty 0.360(3) (± 2 2 \) 1 
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Figure 1: Z^ 07 defined in eq. (31), with w + , w computed at NLO in perturbation theory 
and Z + (fi), Z(fi) computed non-perturbatively in the RI-MOM scheme. The scale \x is 
in GeV. ' 
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The first error is statistical while the second one is the systematic uncertainty estimated 
as explained in sec. 6.2. The conversion to the MS-NDR and MS-HV renormalization 

schemes can be readily done by using eqs. (20)-(23) (we take a^ S (/i = 2 GeV) = 0.1953 
obtained in the quenched approximation at NLO with Aqcd = 0.250 GeV). The non- 
perturbatively obtained results in eqs. (33) and (34) are remarkably similar to the per- 
turbative results in eq. (25) (the similarity is less remarkable if mean field improvement 
is used to estimate the higher order terms in perturbation theory). Moreover, we can 
compare our non-perturbative results, computed with the data sets 1 and 2, with the 
ones obtained previously in ref. [54] at the same f3 but using a smaller volume and data 
set 1 only. In this study, a complete basis of five parity odd AI = 3/2 operators was 
considered. For our purposes we only need to consider three of these, namely (in the 
notation of ref. [51]) Qf, Q + , and <2 + . In fact, with the appropriate choice of flavour 
quantum numbers, there is the following correspondence between the two notations: 

4 = 2g+ O s 7 /2 = 2Q+, 3 8 /2 = -AQ+; 

Z + = Zu, Zj7 = Z 2 2, Z 7 g = —-Z 2 3, Z$8 = Z 3 z, Zg 7 = — 2 Z32, (35) 

and the renormalization structure is therefore the one indicated in eqs. (16) and (17). 
The values obtained in ref. [54] were the following, 

z ri-mom = 0.608(4) (14) 0.604(3) (5) 



rRI-MOM 



0.673(3)(10) -0.139(1)(5) \ / 0.666(3)(3) -0.140(1)(2) \ (36) 
-0.056(2)(2) 0.392(3)(28) I 1 -0.056(2)(2) 0.375(3)(14) J 



where the first corresponds to the small volume 16 3 x 52 and the second to data set 1. 
As can be seen by comparing eqs. (33) and (34) with eq.(36), the agreement is excellent, 
suggesting that finite-volume effects are very small. 



5 The Analysis 



In this section we discuss the analysis of our results, starting with the extraction of the 
matrix elements. We then discuss the determination of the shift of the energy of the 
two-pions as a result of finite- volume effects and finally we explain our procedure for the 
chiral extrapolation. 



5.1 Extraction of the Matrix Elements 



We now explain the procedure we use to determine the matrix elements in the SPQR 
kinematics [14,15], where one of the final-state pions has a non-zero spatial momentum, 
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Figure 2: Plateaus for the extraction of the form factors a 3 .MgpQ R . 



and the other two mesons are at rest. In order to remove the 1=1 component in the 
two-pion state, we take the symmetric combination of the two states |vr + (p)7r°(0)) and 
|7r + (0)7r°(p)), where \p\ = or 2tt/L, i.e. we evaluate the following average: 



MU^M^E^Mk) = - (<7r + (p)7r°(O)|0,|tf + (O)> + <7r+(oV(^l^+(0)>) , (37) 

where E n is the energy of the corresponding pion and Oi are the operators defined in 
eqs. (l)-(3). 

In Euclidean space, up to NLO precision in the chiral expansion and neglecting finite- 
volume effects, the magnitude of -MgpQ R can be extracted from the ratios of correlation 
functions 

1 / C^(t,p) + C^(t,p) 

2 \c K (t 6x , 0)0,(^,0)^,^; 



a 6 ZlZ K ] 



large t 



a 3 M 



(0 

SPQR' 



(38) 



where C+ 1 and are defined in eq.(12) and A^spqr represents the magnitude of 
the matrix element of the bare lattice operator Oi on the right-hand side of eq. (37). 
Zp = (0\(f)p(0)\P) is the wave-function renormalization of the pseudoscalar meson state 
P, and can be obtained from the fit to 



c P (t,o) 



(a 2 Zp) 2 



e -aMp- cogh 



(39) 



at large t, where T is the number of time slices (T=64 in this study). Fig. 2 shows some 
of the plateaus for a 3 A^sPQR- 
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The value of the pseudoscalar decay constant fp corresponding to the (non-degenerate) 
mass Mp is obtained via the one-parameter fit (using Zp and Mp obtained from eq. 
(39)) 



( ^ - t 



(40) 

,Cp P (t,0)J \M P J \Z A 1 

at large t. In order to obtain / we neglect the effects of SU(3) breaking and we extrap- 
olate quadratically on the smallest 5 masses. We obtain / = 0.139(5) GeV from data 
set 1 and / = 0.135(4) GeV from set 2. The renormalization constant Za is determined 
non-perturbatively, as discussed in sec. 4. 

Eq. (38) corresponds to the case in which NLO is applicable and in which 

finite-volume effects are neglected. More generally we use the finite-volume formalism of 
refs. [25,26] to write: 

Cf(t,i?) + Cl?\t,it) \ 3 2 , large t 

la Z_Zk, 



2 \cv(tfix,o)cv(tfix,o)a(t,$ 

{a 3 |^« QR | cosS^W) + O (£)} e-^x, (41) 

where W is the two pion energy in a finite volume so that AW = W — E n — M n and 
5 I=2 (W) is the 1 = 2 two-pion phase shift. We now discuss the finite-volume effects in 
more detail. 



5.2 Two Pion Energy Shift in a Finite Volume 

In this section we discuss the extraction of the energy shift of the two-pion state with the 
SPQR kinematics, AW = W — E v — M n , using a data set of 750 configurations (of which 
set 1 is the subset on which the four-point functions were also computed). We consider 
pions of three masses, each composed of a degenerate quark- ant iquark pair. The masses 
of the pions in lattice units are M n = 0.4438(7), 0.3590(8) and 0.2557(13) corresponding 
to hopping parameters k = 0.13300, 0.13376 and 0.13449 respectively. 

Consider the correlation function 

G4*(t,p) = i((P. + (t,p)P.o(t,O)0i + (O,O)0io(O,O)) 

+ (P^+ (t, 6)P n o (t, p)<Pl + (0 , 0)^o (0,0))) (42) 

where the interpolating operators in momentum and position space, P n and 0^ respec- 
tively, are defined in eq. (13). We have symmetrized over the two interpolating operators 
at fixed momentum P n + and P n o in order to obtain a pure 1 = 2 state. The most general 
four-pion correlation function may be represented in terms of the four diagrams in fig. 3. 
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Moreover, in the SU(2) symmetric limit, since the two interpolating operators at t = 
are at the same spatial point, the symmetrization performed in eq. (42) is automatic. 
The 1 = 2 correlation function in eq. (42) is given by the combination D — C of dia- 
grams in fig. 3, while the corresponding 1 = correlator is given by the combination 
D + C/2 — 3R + 3/2 V. Here we only consider the 1 = 2 channel, which is the relevant 
one for AI = 3/2 transitions. 



To determine the two-pion energy shift AW, we compute the ratio 

G^(t;p) 



R An (t;p) 



C„(t,p)C w {t,Q 



(43) 



where C n is defined in eq. (14), and, in order to isolate the lightest states, we consider 
the behaviour of the correlation functions in the region T ^> t ^> 0. 

On a lattice with periodic boundary condition in the temporal direction G Av has two 
types of contribution; either a two-pion state propagates between and t (or between t 
and T) or there is a single pion propagating between and t and another between t and 
T. In the region T ^> t ^> the ratio is therefore given by: 



R^(t,p) 



X 1 (e 



-wt 



-W(T-t) )+X 



2 le ' E ^e 



-M„(T-t) 



+ e 



-M w t -£; 7r (T-^)^ 



X ie ~ w 2 cosh (W(t - i)) + X 2 e 



cosh((^-M 7r )(t-f)) 



2 X,e- 



cosh (M n (t - I)) cosh (E n (t - |)) 



(44) 



where the constants -^1,2,3 are independent of time but do depend on the finite volume 
energies of the \ir) and states. E n is the energy of a pion with momentum p. We 
have not displayed the contribution of the excited states. 

In fig. 4(a) we plot R^(t,0) and -R 47r (t,p min ) as a function of t for the lightest 
quark mass (corresponding to k = 0.13449). \p min \ = 2n/L and we average over the 
six equivalent directions. In fig. 4(b) we plot R4n(t,p min ) for all three pion masses. 
In order to determine AW we have to fit the R4 n (t,pYs in a range of t such that 
the first term in the numerator of eq. (44) dominates over the second. This condition 




Figure 3: Feynman diagrams which appear in the computation of the most generic four 
pion correlation function as G An {t,p) in the SU(2) symmetric limit. 
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0.4438(7) 


0.3590(8) 


0.2557(13) 


AW f=0 


0.0054(6)(7) 


0.0062(7)(7) 


0.0066(9)(7) 




0.0086(8)(10) 


0.0109(11)(15) 


0.0152(27)(15) 


AW L L 


0.0064(3) 


0.0072(3) 


0.0083(5) 


AW BG 


0.0052(2) 


0.0060(2) 


0.0071(3) 



Table 2: Numerical results (AWp =0 and AW^ = 2 2L ) and theoretical predictions (AWl and 
AWbg for p = 0, defined in eqs. (45) and (48)). For the numerical results, which are 
given in lattice units, we quote statistical and systematic uncertainties in that order. 



is in addition to the usual requirement that the time intervals are sufficiently large 
for the contributions from excited states to be suppressed. In practice we extract the 
energy shifts by fitting R4, n (t,p) to Aexp(-AWt) in the range t E [14 — 15, 22 — 23] 
(symmetrizing the contributions from the correlation functions over times in the forward 
and backward halves of the lattice), with the exact choice of interval depending on the 
mass and momentum. In order to estimate the systematic error, we vary the limits of 
the time interval used for the fit and look at the spread of the values thus obtained. The 
results are reported in tab. 2 (preliminary results were presented in ref. [32]). 

In tab. 2 we also present theoretical predictions for the energy shift obtained using 
the Luscher formula and tree-level in full QCD (AWl) [59-62] and by using one- 
loop quenched x?T (AWbg) [30]. In both cases the predictions are for p = 0. In spite 
of reservations that our pion masses may be too large for xPT to be applicable, the 
agreement of our lattice results with the theoretical predictions from one-loop quenched 
%PT is very good and reassuring. On the other hand, we observe discrepancies between 
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Figure 4: (a) R in (t,p) for aM^ = 0.257 and for both momenta, on the whole lattice, (b) 
R^nityft mm) f° r the three masses of the pion. 
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our results and the predictions of Luscher's formula (which however is obtained in full 
QCD). We now outline the determination of AWl and AWbg- 



The Luscher formula for the energy shift gives 

AW L = W - 2M. = -i^L [i + c q + C2 4| + 0(L - 6) } (45) 



where c x = -2.837297, c 2 = 6.375183 (46) 
and do is the scattering length in the s-wave, 1 = 2 channel. At tree level in full 

°° = ' (47) 

We use the values of M n obtained in our simulation and the corresponding values of 

AWbg is obtained by using one-loop quenched ChPT: 
AW BG = AW tiec + AW 1 ^ = — ^— + M n [B 2 (M n L)5 2 + A 2 (M n L)Se + 0' 



2 



2 f2 L 3 ' ' , -V(M 7r L) 3 - / 

(48) 

where the first term corresponds to the leading term in the Luscher formula after sub- 
stitution of the expression for the scattering length in eq. (47), while the second term 
comes from one- loop diagrams. In the region of large M n L, A2{M n L) and B^iM^V) can 
be approximated by 



Finally 



4s ^S and -iS' (50) 

where m is the 7/ mass. Taking the full QCD estimate vn^j?> = (500 MeV) 2 would 
give 5 ~ 0.18. We vary 5 between 0.08 and 0.18 (which is approximately the range of 
values typically obtained in quenched simulations). In this interval, however, the one- 
loop contribution is small for the parameters of our simulation, and we observe a very 
mild dependence of AWbg 011 $ which is included in the quoted uncertainty. In tab. 2 
we present the results obtained for 5 = 0.10. 

In our data set 2 for K — > tttt decays we have not evaluated the correlation function 
G^ n and therefore cannot determine AW at the corresponding pion masses directly. 
However, as shown in tab. 1 the pion masses in this data set are very close to those 
in set 1 and we therefore obtain the energy shifts at these masses by interpolating or 
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extrapolating the results obtained in set 1. We find that for the pion masses aM w = 
0.4642, 0.3663 and 0.2716, the corresponding aAW p=0 = 0.0052, 0.0061, 0.0065 and 
aAW^s* = 0.0080, 0.0107, 0.0146. 

Finally we use the extracted values of AWp for SPQR kinematics to correct the ratio 
on the left-hand side of eq. (41) for the finite- volume effects e~ Awt ^. Results for the 
matrix elements of bare operators obtained in this way with cos(5) = 1 are reported in 
tab. 3. They still contain finite-volume effects, and in order to remove these one would 
need to know the Lellouch-Luscher factors for SPQR kinematics. These are currently 
unknown and we estimate the corresponding errors in sec. 6.3. 

5.3 Chiral Extrapolations and the Numerical Results 

We start by fitting our lattice data to the formulae of the NLO chiral expansion, including 
the chiral logarithms obtained in ref. [15], determine the low energy constants and use 
these to extrapolate to the physical point. When applicable, we use the numerical values 
for the couplings a = 0.1 and mo = 0.5 GeV in the terms a(d fl &o)(d fl Qo) — vra^o °f 
the quenched chiral Lagrangian, where $o is the super-?]' field. The results are not very 
sensitive to the variation of a and m within a reasonable range. The fits are poor, 
which may not be surprising since we are trying to use the chiral expansion in a range of 
masses and energies which are rather large ([0.52, 1.5] GeV). To improve the situation 
we have also performed fits in restricted intervals, introducing an upper cut-off for the 
masses and energies (we vary this cut-off down to a value of 0.8 GeV). Since the validity 
of quenched %PT (q%PT) may be questioned even in the restricted interval we have also 
performed fits using only the polynomial terms of q%PT (i.e. the NLO expressions from 
q%PT but without the logarithms) which have the same functional form as those in %PT 
of full QCD. 

Of the fits which we have tried the best results have been obtained using the polyno- 
mial fits. This is illustrated in fig. 5, where we display the ratio of the lattice amplitude 
and the fitted amplitude for q%PT and polynomial fits. This exercise shows that the 
NLO polynomial correction are essential to describe the data in the range we have simu- 
lated, whereas the addition of the chiral logarithms (which have fixed coefficients) spoils 
the quality of the fits. Since at sufficiently low masses the NLO chiral logarithms must 
be present, we have also tried the centaur procedure in which we smoothly match the 
polynomial fits (the horse) to the NLO %PT formula (the man) at some kinematical 
point. The details of the matching procedure will be described below. The centaur 
procedure allows us to estimate the systematic uncertainty in the extrapolation of our 
data to the physical point. Of course, these uncertainties will not be fully under control 
until simulations in full QCD are performed at sufficiently small masses so that the NLO 
%PT formula is seen to hold. 

We perform the fits on the two subsets of our data corresponding to the two indepen- 
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Table 3: Numerical values of the bare amplitudes A^spqr(p) = -^spQRl^r, E n {p), Mk) 
as obtained from the ratio on the l.h.s. of eq. (41) multiplied by the FV correction 
e Am fix. M n assumes the values la, lb, lc and Mk respectively {la, 2a, . . .}, {16, 2b, . . .}, 
{lc,2c,...}. 
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Figure 5: Plots illustrating the relative quality of the q;\PT and polynomial fits for (Og). 
The 18 points for each of the two values of momentum correspond to the masses from 
data set 2 in tab. 1 and are plotted as a function of Mk- 

dent simulations: 

1) in set 1 masses and energies have values in the interval [0.52, 1.5] GeV; 

2) in set 2 their values are restricted to be in the interval [0.54, 1.2] GeV. This set 
was generated in order to increase the number of points at smaller masses. 

We start with a discussion of the evaluation of j = 2(7m\O 7j8 \K ) for which we need 
only a relatively small number of low energy constants (five, one at LO and four more 
at NLO). In addition the chiral expansion starts at 0(M® K ) (rather then 0(M% K )) and 
so the fits are much more stable and an estimate of the matrix elements is possible to a 
reasonable accuracy. As explained above, in addition to using NLO qxPT formulae, we 
also fit the data to the polynomial obtained by removing the chiral logarithms. 

The results obtained using the procedures described above are presented in tab. 4. 
(^7,8)/=^ represents the value of the matrix element at physical values of the masses, 
energies and momenta at NLO in the chiral expansion. By (Ot^)1=2 we denote the 
contribution of the leading order to the physical result. Note that this is the value of 
the matrix element in the chiral limit, which is the quantity determined from previous 
studies, in which only K — > it matrix elements were calculated. We observe from tab. 4 
(see also fig. 5) that the polynomial form gives by far the best description of our lattice 
data for / =2 {7T7T | C7.8 1 JC ) , indicating that the data are not in the kinematic range where 
the NLO chiral expansion for these amplitudes are valid. 

From our study we see that the use of qxPT does not describe our lattice data and 
therefore using fits based on q%PT alone can lead to unreasonable results. However, in a 
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fit ansatz 


points 


<0 7 >gf (o 7 )U ttf 


(Os)fj 2 s (Os)} =2 jfa 


Set 1 (340 configurations) 


poly., M K ,E n < 1.5 
poly., M K ,E n < 1.2 
poly, M K , E w < 1.0 


18 
13 
10 


0.130(15) 0.155(16) 0.11 
0.135(15) 0.161(17) 0.13 
0.140(16) 0.161(17) 0.10 


0.746(91) 0.802(99) 0.024 
0.750(93) 0.792(98) 0.020 
0.753(93) 0.790(98) 0.015 


q X PT, M K ,E n < 1.5 
q X PT, < 1.2 
q X PT, Mjr, -E^ < 1.0 


18 
13 
10 


0.203(21) 0.308(32) 14.8 
0.302(32) 0.492(51) 5.83 
0.276(29) 0.477(50) 4.86 


1.37(17) 2.06(27) 8.95 
1.51(19) 2.32(30) 3.46 
1.37(18) 2.20(29) 2.88 


Set 2 (480 configurations) 


poly, M K ,E n < 1.2 
poly, Mtf,^ < 1.0 
poly, M K ,E n < 0.8 


36 
24 
16 


0.127(11) 0.143(12) 0.022 
0.128(12) 0.143(13) 0.014 
0.132(13) 0.142(13) 0.013 


0.734(64) 0.761(70) 0.022 
0.732(68) 0.754(74) 0.019 
0.732(77) 0.753(77) 0.025 


q X PT, Mff, E w < 1.2 
q X PT, Mif,^ < 1.0 
q X PT, Mr-, -Ett < 0.8 


36 
24 
16 


0.268(26) 0.425(42) 13.4 
0.265(28) 0.442(47) 7.71 
0.224(25) 0.412(47) 6.06 


1.53(16) 2.33(26) 11.0 
1.47(16) 2.34(27) 6.36 
1.22(14) 2.18(26) 5.01 



Table 4: Results obtained using various fitting procedures to our data for the matrix 
elements i = 2(Tnr\Oj\K°) and 7=2 { 7T7T | | i^ ) renormalised in the NDR-MS scheme at 
2 GeV. Energies and masses are given in GeV while the matrix elements are in GeV 3 . 
When fitting using the qxPT formula we take a = 0.1 and mo = 0.5 GeV. Notice that 
we perform uncorrelated fits, which is the reason for the very small x 2 /d.o.f. for good 
fits. 

more general context, it has been argued that one may combine the behaviour expected 
in chiral perturbation theory at low masses with the observed polynomial dependence on 
the kinematical parameters in lattice simulations [63,64]. We now attempt to perform 
such a matching. 

For the most general kinematics the amplitudes .M^ 7 ' 8 -* are function of the Lorentz 
invariants p\ + = p£ = M%, p 2 K = M|, p K • Pn+, Pk • Pn° and p n + ■ p n o, where p K , p n + 
and p n o are four- momenta of the kaon and the two pions. The 1 = 2 |7r + 7r°) final state is 
symmetric and therefore the dependence on px ■ p n + and on px • Pn is the same. In the 
region of small masses and energies the matrix element is given by %PT at NLO, which 
for 8 , for example, takes the form: 

•^general 1 = A (x) ( X + [ chiral ^°S] gen (M R , M% , p^+ -pK,Pn° ' PK,Pir+ ' JV)) 

+B 1(x) M K + B 2 ( x )(p 7T + + p n o) ■ p K + B z{x) Ml + B 4 &)(p,r+ ■ Pn°), (51) 
where A( x ) and the B^'s are unknown constants. On the other hand, in the region 
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accessible to our simulation the data are well described by a polynomial of the form: 



•Mgenertf = A (p) + B l(p) M K + B 2(p)(Pn+ + Pifi) " PK + B 3{p )M* + B 4(p) (p n + ■ p w o) . (52) 

Our procedure is to fit the lattice data to eq. (52), extract and the B^s and then 
to impose the following smoothness condition 



M%*£) = (M 

general j gpQR ^ 



(8),poly\ 



) 

m.p. 



general y spQR m p ^ g eneral / SPQR 

(53) 



dM (8) ' xPT \ fdM {8) ' poly 

J 'general 1 | J general 



d0Ji } \ duji 

I SPQR m.p. \ / SPQR m.p 

where the Ui, i — 1 -4, are the four variables M|-, M%, p w + -p^o and Pk'P-k+- The subscript 
"SPQR m.p." indicates that the matching point is chosen in the SPQR kinematics and 
is thus defined by the values of the masses of the pion and kaon and the energy of the 
moving pion (M*,M^,E*) at which the matching is performed. This set of conditions 
determines the values of and -Bj( x ), which can then be used to compute the physical 
amplitude using 

•Mp 8 hys PT = Ax) + [ Chiral l OgUys(m 2 K,ml)) 

+ (b Hx) + B 2{x) + -£ 4(x) ) rn 2 K + (s 3W - B 4(x) ) ml (54) 

Tab. 5 shows the results for / = 2(7T7r|(97 i 8|i^ ) from this centaur procedure, with different 
choices of matching points. An example of this procedure is represented graphically in 
fig. 6. 

In the ratio A4^/Ai^ chiral logarithms and finite- volume corrections cancel since 
the two operators belong to the same representation [65]. This ratio can also be used 
to compare lattice results to other theoretical approaches [66-70]. Results are reported 
in tab. 6. In contrast to the case of (Oj) and ((9 8 ) individually, this ratio is not so well 
described by the polynomial form coming from x?T. This is clear by comparing the 
values for the xVd-o.f of tab. 6 with those of tab. 4 and also by comparing fig. 7 with 
fig. 5. We will discuss these results in sec. 7. 

Implicit in our work is the assumption that, in the region in which we have data, 
the lattice results are reliable even though the simulations have been performed in the 
quenched approximation. When performing the centaur matching we join the polynomial 
fit of our data with the NLO formula corresponding to full QCD. At low masses 
and energies (where we do not have data) the behaviour predicted by quenched is 
very different from that in full QCD. However, in spite of the fact that we have data 
from a quenched simulation it does not make sense to extrapolate to the physical point 
using qxPT. Our central value for / =2 (vr7r|(9 7i8 | J^°) is obtained by the centaur matching, 
with full chiral logarithms, at the matching point (M*, E*, M£) = (0.4,0.4,0.5) GeV. 
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At the same time we also obtain the estimate of the LECs in full QCD reported in tab. 7. 
Other results reported in these tables are used to estimate systematic errors, and will be 
discussed in detail in sec. 6.1. 

We now turn to the matrix elements of the operator O4. There is one LEC at leading 
order in xPT and six more at NLO, so that the fits to the NLO chiral expansion have 
seven parameters. Again we find that our data is not well represented by quenched chiral 
perturbation theory at NLO. If instead we remove the chiral logarithms from the fit, we 
do obtain fits with reasonable values of \ 2 ■ The coefficients of the polynomial terms 
at NLO (which would correspond to the NLO LECs if the chiral logarithms had been 
included) are very poorly determined however, making any extrapolation to the chiral 
limit impossible. On the other hand, the coefficient of the leading term is found to be 
stable, and is significantly lower (by about 40%) than the value obtained keeping only 
LO chiral perturbation theory. Thus the contribution of the NLO terms is clearly visible 
even if, due to the high number of free parameters in the fit, it is not possible to evaluate 
the corresponding LECs with sufficient precision. We therefore do not discuss further 
the evaluation of the matrix elements of O4 beyond showing one indicative plot. In fig. 8 
we plot the values of the bare matrix element -MgpQ R (0) from dataset 2 at the three 
values of the pion mass with degenerate quarks. Also shown are linear and quadratic fits 
to these three points, from which it can be seen that the extrapolated value in the chiral 
limit is close to zero as expected from chiral perturbation theory (indeed the quadratic 
fit gives a value consistent with zero). However, until we have lattice data in the chiral 
regime, plots such as that in fig. 8 are at best qualitative. 



6 Systematic uncertainties 

We now discuss the main systematic errors which affect the values for the renormalized 
matrix elements of the electroweak penguins obtained in this work. 

6.1 Uncertainties from the chiral extrapolation 

Our results indicate that we are not in a regime where qxPT is valid. Completely reliable 
results can only be obtained by performing full QCD simulations in a regime where 
%PT describes the data and can therefore be used for the extrapolation to the physical 
point. In the meantime we use the centaur procedure of matching the polynomials which 
describe our data to the NLO chiral expansion for (7r + 7r° | (f^^g 2 • This procedure is 
valid under the assumption that there is an overlap region in which both descriptions 
hold. We estimate the corresponding systematic error by varying the matching point 
below our data down to the chiral limit (the results are insensitive to varying the cut-off 
on the masses and energies of our data points). 
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0.779(68 
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Table 5: Results for i = 2{nir\O'i i s\K ) via xPT-polynomial matching performed with full 
and quenched chiral logarithms from Set 1 (340 configurations), first block of results, and 
Set 2 (480 configurations), second block. The matrix elements are renormalised in the 
NDR-MS scheme at 2 GeV. Energies and masses are in GeV while the matrix elements 
are in GeV 3 . In the q%PT formula we use a = 0.1 and m = 0.5 GeV. 
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Figure 6: xPT-polynomial matching for (0 8 ) I=2 (in GeV 3 ) at (M*,E*,M^) = 
(0.4,0.4,0.41) GeV (we plot only data points with M n = E n = M K ). For reasons of 
numerical stability, it is convenient to use a matching point in which is offset from 
M* and M*. The point labelled as chiral limit corresponds to M n = E n = Mk = 
and that labelled as physical point corresponds to the physical values of the masses and 
energies, and hence does not lie on the SPQR curve. 



data sets 


points 


mfsimri, si? 


Set 1, M K ,E n < 1.5 
Set 1, M K ,E V < 1.0 


18 
10 


0.157(12) 
0.171(11) 


H 


3.11 
3.02 


Set 2, M K ,Ev < 1.2 
Set 2, M K ,E W < 0.8 


36 
16 


0.153(10) 
0.163(10) 


r< 
H 


1.33 
1.14 



Table 6: Results of the polynomial fit for the ratio / Ai^ where the operators are 
renormalised in the NDR-MS scheme at 2 GeV. 
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Figure 7: Quality of the polynomial fit for (Or)/ (0$). 



LECs 


Set 1 O r 


Set 2 C> 7 


Set 1 O s 


Set 2 8 


7// 3 [GeV 3 ] 


0.0705(75) 


0.0653(54) 


0.365(45) 


0.347(32) 


{ATTff(8 1 +5 2 )h 


-0.42(90) 


-1.08(59) 


-1.20(64) 


-1.64(53) 


(4vr/) 2 (|-| + 5 3 )/7 


-0.46(33) 


-0.36(31) 


-1.27(20) 


-1.35(19) 


(4vr/) 2 (54 + 5 5 )77 


0.12(43) 


-0.06(24) 


0.66(26) 


0.44(21) 




-0.59(18) 


-0.39(13) 


-0.92(13) 


-0.75(11) 



Table 7: Results for the LECs of i = 2(7m\Oi t s\K°) (where operators are renormalised 
in the NDR-MS scheme at 2 GeV) from the centaur matching with (M*, E*, M£) = 
(0.4, 0.4, 0.5) GeV. Only statistical errors are displaied. The values for / are 
0.139(5) GeV from set 1 and 0.135(4) GeV from set 2. For the definitions of (5 4 + S 5 )' 
and Sq see eq. (11). 
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Figure 8: Chiral behaviour of the bare amplitude (in lattice units) A^spqr(O) f° r pions 
with degenerate quarks (lines la, lb, lc of Tab. 3 for dataset 2) versus M%. The numbers 
in the box correspond to the extrapolated values in the chiral limit. 
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6.2 Uncertainties from operator matching 

The main source of systematic error in the evaluation of the renormalization constants is 
the residual dependence on the renormalization scale \i of Z+ GI and Z in eq. (31). As 
explained in sec. 4.3, this residual dependence is probably due to higher order corrections 
not taken into account by the NLO perturbative evaluation of the evolution functions w + 
and w. Since the scale /i = 2 GeV at which we evaluate the renormalization constants 
lies inside the range of momenta at which we compute the renormalization constants, we 
estimate the systematic error in the following way: we take the renormalization constants 
computed non-perturbatively at a given scale /i and we run them up or down to the 
required scale fi. The systematic error is then estimated from the deviations of the 
results obtained by varying fiQ in the whole range in which we compute numerically the 
renormalization constants (i.e. fio G [1.6, 2.75] GeV). In order to reduce this uncertainty 
we would need to perform simulations at larger values of a" 1 . 

6.3 Uncertainties from the 7r— 7r phase shift and the finite- 
volume effects 

In this study we have used matrix elements computed with the SPQR kinematics for 
which the two-pion system is not at rest. The Lellouch-Liischer (LL) factor relating 
the finite-volume matrix elements to physical amplitudes is valid in the centre-of-mass 
frame, and the corresponding formula for the moving frame is currently unknown. This 
introduces a systematic error and underlines the importance of developing a theory of 
finite-volume effects for two-meson states in a moving frame. Using the LL factor as a 
guide with phase-shifts estimated using we expect these corrections to be of order 
10-15%, but their effect on the chiral extrapolation is difficult to estimate. 

As shown in sec. 5.1 and 5.2, we were able to eliminate the systematic error coming 
from the shift AW of the two-pion energy in a finite volume. After applying this correc- 
tion, since we divide the correlator C+ + Cq^ by two non-interacting pion propagators 
(see eq. (38)), we obtain the matrix element multiplied by the cosine of the phase-shift, 
cos(S I=2 (W)). For the data points in which both pions are at rest, 5 I=2 (W) = and 
there is no correction. For those where one of the two pions has non zero momentum 
(i.e. in a moving frame), a first attempt to study the relation between the finite- volume 
energy shift and the infinite volume phase-shift has been presented in ref. [71] but a thor- 
ough understanding of this problem proves to be very involved and might not be resolved 
soon. For this reason we give here only an approximate estimate of cos(5 I=2 (W)) based 
on the computation of the phase shifts obtained with Nf = 2 dynamical simulations in 
the laboratory frame [72] 6 . For all of the pion masses in our simulation cos(5 I=2 (W)) 

6 In view of the fact that the phase shifts for pions in the center of mass frame are very similar in 
the dynamical and in the quenched case at the same lattice spacing and on the same volume [72.73], 
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turns out to be larger than 0.98. 



6.4 Quenching effects 



A reliable estimation of the quenching effects can only be obtained when the unquenched 
(or partially quenched) lattice data are available. The matrix elements presented in this 
paper have dimension three, therefore suffer from the uncertainties in the determination 
of the lattice spacing, which is one effect of the quenched approximation. By performing 
phenomenological studies using ratios of matrix elements it may be hoped that the 
uncertainty is reduced, but it must be remembered that quenching errors remain and 
that simulations with dynamical quarks will be required to eliminate them. 



7 Final results 

Our final results from the two data sets (at the precision of NLO in the chiral expansion) 



arc 



I=2 (nn\0 8 \KXt y ° s = 0.691(84) (tg) GeV 3 

,NLO 
'phys 



/=2 <Ha|^°>pts = 0.119(14) 1I 2 tl 1 ) GeV 3 (55) 



from set 1 and 

/=2 <H(?8|A^ h L ° = 0.682(59) (+S) (tlfj GeV 3 
i^(™\Or\K°)X = 0.117(11) (jj) (i?) GeV 3 (56) 

from set 2, where the matrix elements are renormalized in the NDR-MS scheme at 
2 GeV 7 . The first error is statistical, the second is due to the uncertainties in the 
non-perturbative matching to the continuum renormalization scheme and the third is 
due to the xPT-polynomial matching. In estimating the third error, we have taken the 
central value as that obtained by matching to full chiral logarithms at (M*,M^) = 
(0.4 GeV, 0.5 GeV), the lower limit to be the value obtained by matching at (0.4 GeV, 
0.41 GeV) and the upper limit to be the value obtained from the polynomial fit. With 

we assume that this is the case also in the laboratory frame, for which only dynamical studies exists. 
7 The normalization of states and definition of operators in the Introduction implies that 

/=2 <H0 7 , 8 |A"°) - y|(7r+7r |O 7 , 8 |^ + ) = sJ\(k + AO^\K+). 
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this choice, the error bar spans the range of values in the top half of table 5. The results 
from the two data sets are consistent and as our best estimate we take 

I=2 (nir\0 8 \KXky s = 0.68 (6) (4) (5) GeV 3 

/=2 <H^°>?h° = 0.12 (1) (1) (1) GeV 3 . (57) 

Taking the ratio of these two matrix elements directly at the physical point we obtain 
the values 

0.173(17) (±1) (±1) from set 1 

0.172(14) (tfj (t{) from set 2. (58) 

The last error is again due to the uncertainty in the xPT-polynomial matching. From 
the last column of tab. 5 one can appreciate how this ratio is independent of the effects 
of the chiral logarithms, i.e. from the choice of the matching point (M*,E*,M^) and 
from the use of xPT or qxPT formulae. The results in eq. (58) can be compared with 
those obtained from the polynomial fit of the ratio A4^ / 'A4^ in tab. 6 (which has a 
relatively poor x 2 )- By reducing the range of masses and energies used in the fit, its 
quality increases and the values for the ratio get closer to those reported in eq. (58). 

The contribution of the leading order term to the matrix elements, which also corre- 
sponds to their value in the chiral limit, is 

I=2 (nn\O 8 \K ), = 0.843(104) (+f 7 ) GeV 3 

I=2 (7rn\0 7 \K% = 0.163(17) (±» s ) (+ 9 7 ) GeV 3 (59) 

from set 1 and 

I=2 (nn\0 8 \K°), = 0.800(74) (^) GeV 3 

I=2 (7i7i\O 7 \K ), = 0.151(13) GeV 3 (60) 

from set 2, where the sources of errors are the same as in eq. (55). From eqs. (55) and 
(59), it is clear that the NLO contribution of the chiral expansion is significant. It tends 
to decrease these matrix elements by about 15 — 25%. In the kinematical range of the 
lattice data these NLO contributions are very significant. In fact, had we computed the 
matrix elements at LO in %PT (i.e. by fitting our data points to a constant) we would 
have obtained the following results: 

I=2 (7m\0 8 \K%° s = 0.951(74) GeV 3 

I=2 (7in\0 7 \K%L = 0-082(9) GeV 3 (61) 



/=2(7T7r 


W| A /phys 


1=2^71 
/=2(VT7T 


Os\K°)X° s 


1=2 (tTVT 


Os\KXt° s 
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(07>j = 2 [GeV 6 \ 


{0 8 y I= 2 [GeV 5 ] 


{07)1=2/ (°s)l=2 


This work 


0.16(3) 


0.82(15) 


0.19(3) 


Donini et al. [74] 


0.11(4) 


0.51(10) 


0.22(9) 


RBC [10] 


0.27(3) 


1.1(2) 


0.25(5) 


CP-PACS [11] 


0.24(3) 


1.0(2) 


0.24(6) 


Friot et al. [70] 


0.12(2) 


2.00(36) 


0.057(15) 


Cirigliano et al. [68] 


0.16(10) 


2.22(67) 


0.07(5) 


Bijnens et al. [67] 


0.24(3) 


1.2(8) 


0.20(14) 


Cirigliano et al. [69] 


0.22(5) 


1.50(27) 


0.15(4) 


Narison [66] 


0.21(5) 


1.40(35) 


0.15(5) 



Table 8: Numerical results for the matrix elements of Oj and 0$ (renormalised in the 
NDR-MS scheme at 2 GeV) in the chiral limit from various methods (lattice, dispersion 
relations, large A^+MHA, sum rules). 

from set 1 and 

I=2 (7nr\0 8 \KX°y S = 0.979(45) GeV 3 

I=2 (Tnr\0 7 \K%° s = 0.082(9) GeV 3 (62) 

from set 2. 

In tab. 8 we compare our results with those obtained using different methods (includ- 
ing lattice simulations, dispersion relations, large iV c +MHA and QCD sum rules). None 
of the lattice estimates includes the quenching error. The three previous lattice com- 
putations [10,11,74] used K — > ir matrix elements and at leading order to obtain 
K — > 7T7T amplitudes in the chiral limit. The present work is the first computation with 
two pions in the final state which allows the use of x?T at NLO to extrapolate to the 
physical point and to the chiral limit. All of the lattice computations have been per- 
formed at fixed value of the lattice spacing, but with different actions, different methods 
of computing the renormalization factors and different method of performing the chiral 
extrapolation, making it very difficult to compare the systematic uncertainties and to 
understand the source of the discrepancy among different results. In the next generation 
of computations it will be important to perform simulations at several values of the lat- 
tice spacing in order to study the approach to the continuum limit, as well as going to 
lighter values of the quark masses (the relatively high- values of masses used in current 
and previous simulations is one of the major limitations). 
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The range of values for (Og)} =2 obtained using lattice simulations is lower than that 
obtained using other techniques (although they do overlap within errors). For (Oi)} =2 
there are discrepancies between the different lattice simulations and similar discrepancies 
using other techniques. Finally, one of the main conclusions of the present study is that 
the NLO contribution of at the physical point decreases the value of the matrix 
elements by a 15 — 25%. This is a physical effect which will affect all the estimates 
presented in tab. 8. 



8 Conclusions 

It can be argued that K tttt decay amplitudes represent the Holy Grail of lattice 
simulations, since so many theoretical and technical difficulties must be overcome before 
they can be calculated reliably. The aim of this study was to get a better understanding 
of the practicability of computing correlation functions with two-pion final states. This is 
necessary for extending calculations beyond leading order in chiral perturbation theory. 
Specifically in this paper we have presented the first attempt to calculate the matrix 
elements for K — > nn decays with an / = 2 final state and to determine all the low 
energy constants which appear at next-to-leading order in the chiral expansion. We 
were able to carry out all the necessary steps but with varying degrees of precision and 
success. Our best estimates for the matrix elements of the electroweak penguin operators 
are presented in sec. 7 (see, for example, eq. (57)). In addition to the quoted errors 
should be added those due to the use of the quenched approximation and finite-volume 
corrections. Although we could also determine the matrix elements of the operator O4 
at the masses which we used in the simulation, we were unable to determine the seven 
low-energy constants with sufficient precision to perform the chiral extrapolation. For 
all three operators we conclude that the NLO terms are significant. 

The positive aspects of our exploratory study include: 

1. For the 1 = 2 channel, the finite- volume energy shift can be determined from the 
linear behaviour of the ratio R^ n defined in eq. (43). 

2. The K — > txtx matrix elements can be determined with good precision at the values 
of the masses used in the simulation. 

3. The non-perturbative renormalization of the operators can be performed success- 
fully. 

In addition to removing the quenched approximation by performing simulations with 
dynamical quarks, we conclude that two major improvements are needed before reliable 
and precise results can be obtained using the techniques of this paper: 
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1. Simulations have to be performed in a region of light-quark masses where ;\RT 
is valid at NLO for K — > nir matrix elements and which is sufficiently large that 
the low energy constants can be determined precisely enough for the extrapolation 
to the physical point to be possible. We have seen that this is not the case for 
the masses used in this simulation and for this reason we had to use the centaur 
procedure for the chiral extrapolation. This needs to be eliminated. 

2. A theory of finite- volume corrections needs to be developed for two-pion states 
with a non-zero (total) momentum. For the spectrum this issue was studied in 
ref. [71]. For alternative procedures to determine the low-energy constants at 
next-to-leading order, which exploit more the freedom to vary the quark masses 
in lattice simulations rather than introducing a non-zero momentum for the irn 
system see ref. [75]. 

Some minor improvements to our study can be made relatively easily. For example, the 
factor of cos(5) due to the two-pion sink (and the corresponding finite-volume effects) 
can be eliminated by computing four-pion correlation functions. This would also provide 
a consistency check on the value of the finite- volume energy shift. 
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A Mass Dependence of K — > tt Matrix Elements. 

In this appendix we consider the chiral behaviour of a simpler quantity than the K — > tttt 
matrix elements described in the body of the paper, viz. the K — > ir matrix elements from 
the quenched simulation of ref. [11]. For illustration we consider the matrix element of 
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the operator O s and show that the uncertainties due to the extrapolation in the quark 
mass are also significant in this case. 

In figs. 9 and 10, we present the results obtained using different ansatze for the chiral 
fits to the quenched data for the ratio {7r + \Og\K + )f 7r /f M (taken from ref. [11]). The 
matrix element (n + \Ol /2 \K + ) is determined for M K = M n = M M . f w « 132 MeV is the 
physical pion decay constant, and ]u is the decay constant of the pseudoscalar meson 
with the mass Mm- In order to estimate the physical result we perform the following 
procedures: 

1. We fit the lattice data to a quadratic function of the form + A < 2* > M M + A < f S> M M . 
The fit is shown in figs 9 and 10. 

2. We also fit the data to the prediction of x?T of the form A^\l + [chiral log]) + 
A^Mlj. The chiral logarithms in full QCD for {-K + \Ol /2 \K + ) can be obtained 
from ref. [76]. In the quenched theory, the corresponding chiral logarithms have 
been calculated in ref. [77]. These fits in full and quenched QCD are also shown 
in figs. 9 and 10 respectively. 

3. Finally we carry out the centaur procedure of matching the polynomial fits to the 
lattice data (which again are good representations of the data) to chiral perturba- 
tion theory at lower masses. In the same way as for the K — > -nix matrix elements 
discussed in sec. 5, we assume that there is a region where both the polynomial and 
%PT behaviour are valid so that by matching the two forms at some scale (which, 
in figs. 9 and 10 is chosen to be 465 MeV) we can obtain the coefficients A[ of 
the chiral expansion from the fitted coefficients A± . We perform the matching by 
equating the value of (ti + \0^ 2 ^\K + ) f n / flj and its first derivative with respect to 
M M at the matching scale. 

These figures show that the choice of the ansatz for performing chiral extrapolation can 
indeed result in large systematic errors, because the lattice simulation is carried out at 
relatively heavy quark masses. See for example fig. 9, where the matrix element in the 
chiral limit using the polynomial fit and extrapolation is about twice that obtained with 
the centaur procedure. 
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Figure 9: K — >• 7r matrix element using data from ref. [11]. The fit corresponds 

to full QCD. The vertical line indicates the value of Mm = 465 MeV where the 
polynomial matching is performed. 
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Figure 10: K — > 7r matrix element using data from ref. [11]. The xPT fit corresponds 
to quenched QCD. The vertical line indicates the value of Mm = 465 MeV where the 
xPT-polynomial matching is performed. 
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